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Abstract —Energy harvester based cognitive radio is a promis¬ 
ing solution to address the shortage of both spectrum and energy. 
Since the spectrum access and power consumption patterns are 
interdependent, and the power value harvested from certain 
environmental sources are spatially correlated, the new power 
dimension could provide additional information to enhance the 
spectrum sensing accuracy. In this paper, the Markovian behavior 
of the primary users is considered, based on which we adopt a 
hidden input Markov model to specify the primary vs. secondary 
dynamics in the system. Accordingly, we propose a 2-D spectrum 
and power (harvested) sensing scheme to improve the primary 
user detection performance, which is also capable of estimating 
the primary transmit power level. Theoretical and simulated 
results demonstrate the effectiveness of the proposed scheme, in 
term of the performance gain achieved by considering the new 
power dimension. To the best of our knowledge, this is the first 
work to jointly consider the spectrum and power dimensions for 
the cognitive primary user detection problem. 

Index Terms —Cognitive Radio, Energy Harvesting, Spectrum 
Sensing, Hidden Markov Model, Power Sensing, 2-D Sensing 


I. Introduction 

With the rapid development of wireless services in the 
past few decades, spectrum resource, which is vital and hard- 
limited, faces a critical situation of being scarce. However, 
studies have revealed that we are wasting the spectrum as most 
allocated bands are being under-utilized. To address this prob¬ 
lem, researchers have proposed the idea of dynamic spectrum 
access, which could help increase spectrum efficiency. 

Cognitive Radio (CR) is the well-accepted technology to 
achieve dynamic spectrum access, with its core idea of al¬ 
lowing Secondary Users (SUs) to access spectrum when the 
licensed Primary Users (PUs) are idle. The goal for CRs is 
to maximize the overall spectrum efficiency while preventing 
harmful interference to PU transmissions. One crucial building 
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block of CR is spectrum sensing, which determines whether 
certain spectrum is occupied by some active PUs. 

Many statistical methods have been adopted for spectrum 
sensing in cognitive radio design. The authors in 111 apply 
the Neyman-Perason lemma to study both the local and 
cooperative PU detection schemes, in which the sufficient 
statistics is compared against a certain threshold to detect the 
channel status. When the PU transmission signaling is known 
at the SU side, cyclostationary features could be explored for 
PU detection 0. For wideband cognitive radio systems, in 
addition to energy detection 0, compressive sensing 0 is 
also adopted to efficiently identify spectrum holes. However, 
most of methods above are sensitive to the shadowing/fading 
effects over the PU-to-SU link. 

On the other hand, there are many related works that explore 
the PU Markovian behavior. The existence of PU Markov 
patterns is validated in 0 by real-time measurements in the 
paging band (928-948 MHz), with different effects of false 
alarm and miss detection are studied in (6), where a modified 
forward-backward detection algorithm is proposed to minimize 
the detection risks. In (7), a Hidden Bivariate Markov Model 
(HBMM) is used to quantify both the channel status and its 
dwelling time. For collaborative spectrum sensing, a parameter 
estimation algorithm with classification method is introduced 
in HI to identify the malicious users based on a Hidden 
Markov Model (HMM). 

Meanwhile, due to the growing demand of energy efficiency, 
the energy harvesting based cognitive radio network emerges 
to both improve channel utilization and meet the requirement 
of green communications. In the power domain, Markov 
chain models have been widely accepted mna to specify 
the energy arrival process, for which some other statistical 
models have also been applied. For example, a Poisson model 
with a known intensity Ao is adopted in m and a Gamma 
distribution model is adopted in lfl2ll . Meanwhile, for wind 
energy harvesting systems, Weibull distribution fl3l is widely 
adopted to forecast the wind speed and the corresponding 
harvested power level. 

Unlike traditional communication systems, for the ones 
powered by the environment energy harvesters, power be¬ 
comes a multiple access medium since the power usages across 
different users (especially geographically neighboring ones) 
are correlated CM3. Accordingly, the spectrum access 
and power access are actually correlated events, and such 
correlations could be explored to help with improving the 
PU detection performance in CR, which used to be solely 
dependent on spectrum sensing. In this work, we consider 
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the PU detection problem in energy harvesting based cogni¬ 
tive radio networks. Given the correlation between spectrum 
and power usages, and by considering the spatial correlation 
among energy harvesting users, we propose a 2-D sensing 
scheme that could infer the PU behaviors by jointly learning 
the spectrum and power access dynamics. This is promising 
to improve the PU detection performance since the traditional 
methods, which are solely based on spectrum sensing, have 
certain limitation: When the PU-to-SU transmission is under 
fading/shadowing, the detection performance degrades sharply 
as the channel observation is no longer reliable. Since the PU- 
to-SU channel quality does not affect the power-dimension 
inference, the proposed 2-D scheme could overcome the effect 
of channel fading/shadowing and provide a more reliable 
combined sensing result. 

In addition, traditional spectrum sensing methods only focus 
on sensing the “on-off' activity of PUs. There is also a 
growing demand of knowing the transmission power levels 
of PUs. For example, in m, a novel method that utilizes 
not only the temporal but also the spatial spectrum holes is 
proposed, which requires the knowledge of PU transmission 
power to estimate the PU coverage area. In HD. an optimal SU 
power allocation scheme is provided, which also depends on 
the knowledge of multiple PU transmission power levels. As 
a side product, the 2-D sensing scheme proposed in this paper 
could sense the spectrum and estimate the PU transmit power 
level simultaneously, which could provide more potentials to 
enhance the performance of the energy harvesting based CR 
networks. 

The rest of this paper is organized as follows. Section II 
describes the basic system model. In Section III, the Hidden 
Input Markov Model (HIMM) is proposed to further specify 
the system dynamics. In Section IV, the proposed 2-D sensing 
scheme is analyzed. In Section V, the parameter learning 
algorithm for HIMM is studied. In Section VI, the simulation 
results are provided to show the advantage of the proposed 
scheme. Finally, Section VII concludes the paper. 

II. System Model 

We consider a simple cognitive system with one PU and 
one SU, which are both powered by harvested energy from the 
same renewable source. If the PU does not transmit at certain 
time slots, the harvested energy is discarded since no battery 
is assumed m. In our setup, it is possible that even when 
the PU has data, the energy level may not meet the reliable 
transmission minimum requirement, which implies that the PU 
does not occupy the channel. This is the extra PU idle case 
compared with the traditional PU system where PU is idle only 
if it has no data to transmit. Specifically, let Hq denote the 
case when the PU does not transmit and Hi denote the case 
when the PU occupies the channel; then a formal definition 
for Hq and H\ in an energy harvesting based cognitive radio 
network is: 

Ho : No Enough Energy or No Data Available ; 

Hi : Enough Energy and Data Available. 

Let {E t ,t = 1, 2, 3, • • • } denote the energy arrival process 


at the PU on a time-slotted basis, where each E t could take 
value from a set with a finite cardinality = {i,i + l,i + 
2, • • • ,i + L — 1}, i > 0, i.e., the harvested energy value is 
quantized into L levels. Here L is set as a relatively large num¬ 
ber, such that the effect of discretizing the PU energy level is 
negligible. The channel occupancy state {Ct, t = 1,2,3, • • ■ } 
is a discrete-time process such that each state takes value in 

= {0,1, 2, • • • , M — 1}, with M = 2 in this paper, which 
implies: When Ct = 0 the channel is idle, and when Ct = 1 
the channel is occupied by the PU. 

Let E^ be the minimum energy level required for reliable 
primary transmission and Pq be the probability of no data 
available. Then the relationship between the channel state 
{C t } and the PU energy level {E t } could be expressed as 

P(C t = 0) = P(E t < Eh) + P(E t > E h ) ■ P 0 , 

P{C t = 1) = P(E t > Eh) ■ (1 - Po)- (1) 

Note that we assume no battery installed such that when 
a PU transmits, it uses up all the harvested energy available 
at that time slot; when a PU does not transmit, it discards 
the harvested energy in that time slot. On the SU side, 
since the channel state and PU energy level are not directly 
observable, SU could only estimate the hidden states with its 
available observations. First, as both PU and SU are powered 
by harvested energy from the same renewable source, the 
harvested SU energy could be treated as an observation for the 
PU energy level, due to the spatial correlation of the energy 
harvesting processes. The relationship between the latent PU 
harvested energy E t and the SU harvested energy U t will be 
further discussed in Section III (See Fig. |TJ). 

Assuming that the PU and SU operate on a synchronous 
time-slotted fashion, the sampled received signal from PU to 
SU at time slot t is given as (under real-valued signaling 
assumption) 

H 0 : x t {n) = u t {n ), n = 1,2, ■ ■ • , IV, 

Hi: x t (n) = h t ■ s t (n) + u t (n) , n=l,2, ■■■,N, 

( 2 ) 

where n is the sampling index with a total of N samples 
per slot, St{n) is the signal transmitted by the PU, h t is the 
channel gain constant over each slot, and ut(n) denotes the 
i.i.d. Gaussian noise. The signal sample st(n) is assumed i.i.d 
and independent from the noise. Denote Y t as the summation 
of N received signal energy samples at the SU, i.e., 

JV-l 

Y t = x t (*)■ (3 > 

i=0 

In the following sections, we will first model the mathemat¬ 
ical structure of the SU observation signals and then learn the 
parameters of HIMM, based on which the PU hidden states 
are estimated. 

III. Hidden Input Markov Model 

Markov chain has been widely adopted 00! to specify 
the environmental energy harvesting process, and to model 
the PU data arrival process in traditional cognitive radio 
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Fig. 1. Hidden Input Markov Model Structure 


networks 00. In our work here, we first adopt two discrete¬ 
time Markov chains to represent the PU data arrival and 
energy harvesting processes respectively, then quantify the 2-D 
signaling structure with a HIMM as shown in Fig. [I] 

From Fig. |T] we see that there are several differences 
between the traditional HMMs and the Markov structure 
abstracted in our problem. If the PU energy level is known, 
which implies that Ut reflects E t perfectly, then the structure 
of our problem becomes similar to the Input-Output Hidden 
Markov Model (IOHMM) fl9l . For IOHMM, the training 
process is a supervised learning problem, such that with the 
training input and output, the functional mapping between the 
input and the output could be inferred. Our structure is more 
complicated than IOHMM, since the unobservable Markov 
state Ct is a function of both the previous state Ct-i and the 
input E t , while the input E t is also a hidden Markov process 
from the SU’s point of view. On the other hand, the hidden 
input could be observed not only from the SU energy level Ut, 
but also from the channel observation Y t when the channel is 
active. We thus call this structure as HIMM, for which the 
existing HMM results could not be directly applied. 

As both {E t } and { C ' t } are first-order Markov chains, 
according to (jTJ we have C t = f{C t -i,E t ). Let C l denote 
a collection of channel states from time 1 up to time t (with 
a similar definition for E f ). The first-order Markov property 
implies the following relationship 

P(C t \C t - 1 ,E t )=P(C t \Ct- 1 ,E t ), 

P{E t \E t ~ 1 ) = P{E t \E t -{). (4) 

Throughout the paper the notation P{A) always refers to the 
probability of the random variable taking a particular value, 
i.e., P(A) = P(A = a), unless specified otherwise. 

We use A and B to denote the transition matrices for the 
processes of {E t } and {C t }, respectively. Clearly, as {E t } 
is a traditional Markov chain and takes values from L states, 
A = {a,ij} is a L x L matrix, in which each component a,j 


is given as 

dij = P{E t = j\E t -i = i ), i,j G ££. (5) 

The transition matrix of B = {B(g), q G if} is actually a 
set of matrices. For each q, B (q) is a M x M matrix. Inside 
B(g), each component h tj q indicates the following transition 
relationship of channel states: 

bij, q = P{C t =j\C t -i =i,E t = q), i,j G^,qG^. 

( 6 ) 

Besides the transition probability, the initial probability 
distribution is also important in describing a Markov chain. 
Here we use vector n E and matrix tt c to specify the initial 
distributions, in which each element stands for 

7rf = P(Ei = i), i G 2zf 

trg = P(C! = j|£i =i), i G Sf,j G (7) 

For the HIMM shown in Fig. |T] the SU observations include 
the SU energy level U t and the SU received signal energy Y t . 
Mathematically, due to spatial correlation, we could model Ut 
as a function of E t . Without loss of generality, we assume 
U t G -Sf. Let the L x L matrix D = {dij} be the emission 
matrix with each defined as 

d ij = P(U t =j\E t = i), i,j G . (8) 

Note that conditioned on E t , the only randomness left at Ut is 
the measurement noise, which we could assume independent 
over time. Therefore, we could have the following decompo¬ 
sition: 

t 

P(U t \E t )=l[P(U i \E i ). (9) 

i=l 

On the other hand, based on the Central Limit Theorem 
(CLT), when the number of received signal samples N is 
relatively large, the sum in (|3]i follows a Gaussian distribution. 
Therefore, conditioned on the hidden states, we assume that 
the observation {Y t } follows a Gaussian distribution, having 
its relationship with E t and Ct shown in Fig. |T| According to 
the SU received signal in (|2), when the channel state is idle, 
the output Y t only reflects the energy of channel noise, which 
contains no information about the PU energy level; when the 
channel is busy, the channel output includes the PU signal, 
such that the distribution of Y t is affected by both the hidden 
states Et and C t . For each t, Y t is distributed as 

P(Y t \C t = i,E t = j ) ~ N^al), i G G Se. (10) 

Based on our previous discussion, when * = 0, Hij a 'id erf 
are respectively identical for all j’s. Let matrix fi = {} 
and matrix er 2 = {erf}, the emission probability of Y t is 
then specified by fi and er 2 . Also, given the hidden states, the 
channel observation Y t is conditionally independent over time. 

The initial probability, transition probability, and emis¬ 
sion probability are three important parts that constitute the 
parameters of this unobservable Markov chain. Let 77 = 
{ 7 t e , 7V C , A, B, D, /X, er 2 } be the collection of these param¬ 
eters. When rj is decided, the entire structure of this Markov 
model is specified. Therefore, learning these parameters is 
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Fig. 2. Structure of the Proposed Algorithm 

crucial in order for us to explore the application of this 
unobservable Markov chain, which is used in our PU detection 
scheme. 

The structure of the overall scheme is provided in Fig. [2] to 
better illustrate the work. In the following two sections, we 
first discuss the PU detection scheme assuming an estimated 
parameter ?/; then we present a parameter learning algorithm. 
The arrangement will make the flow more smooth since the 
parameter learning algorithm is quite lengthy. 

IV. 2-D Sensing 

In this section, we analyze the PU detection scheme with 
2-D sensing over the HIMM model. Let t denote the current 
time slot; since HIMM takes the past observations into consid¬ 
eration, at time t the past and current observations { U t , Y f } 
are all available to contribute to the detection process. With 
the estimated parameter vector rf, the probability measure of 
the current PU hidden states could be obtained by recursively 
marginalizing all the past states, while the observations in each 
time slot work as correctors to modify the previously predicted 
result. 


A. Sensing Algorithm 

First, we state the 2-D sensing scheme as follows. The 
detection for the hidden PU channel state C t and the estimation 
for the PU energy level E t could be jointly achieved by 
maximizing the conditional probability of the current hidden 
states, i.e., ( Ct,E t ) is decided as 

C t ,E t = argmax P(C t = c t , E t = e t |f7‘, Y f ). (11) 

c t ,e t 

In CR, such an estimation result further decides the activity 
of SU: When Ct = 0, SU considers the channel as inactive 
and utilizes it; when C t = 1, SU considers the channel as 
occupied and stays silent or transmits carefully based on the 
knowledge of E t . 

In fact, the joint probability P(U t ,Y t ) is the proportional 
constant between P(C t , EtlU 1 ,Y t ) and P(C t , E t , U*, Y f ), 


i.e., by Bayes rules there is a relationship 

P(Ct, Et\U\ Y 4 ) = P(C t , E t , U\ y*)/P(CA Y*) 

c xPiCuEuU^Y*). (12) 

This result implies that the estimated hidden states Ct, E t also 
maximize the joint probability P{Ct, E^U*,Y t ). In ff20l . it 
is demonstrated that such estimated results are based on the 
Maximum A Posteriori (MAP) principle, which is considered 
as an optimal estimate to minimize the Bayes risk for a “hit- 
or-miss” cost function. Next we discuss how to evaluate and 
solve <UD- 

We first give the _ procedure to recursively calculate 

P(C t , EtlU*,Y t ) in ( fill , which is divided into the follow¬ 
ing three steps. Since the Markovian property could not 
be easily applied to decompose the conditional probability 
of the hidden states for recursive calculations, in the first 
step, the probability of the current observation given the past 
P{Ut,Y t \U t ~ 1 ,Y t_1 ) is used to multiply the target condi¬ 
tional probability P{Ct, E^U 1 ,Y t ), such that the production 
after multiplication could be directly decomposed and com¬ 
puted as 

p(u t , Y t \u y t “ 1 )P(c t , Et\u\ y‘) 

= Y Y PiC^Ct-!, Et^t^UuYtlU 1 - 1 ,^- 1 ) 

Ct-1 E t _ i 

= EE P(Y t \Ct,E t )P(U t \E t ) 

Ct-i E t -1 

= P(Y t \Ct,E t )P(Ut\E t )- 

■ v 1 

corrector 

Y Y P{Ct\Ct-i,Et)P(Et\E t - 1 )P{Ct-i,Et-i\U*- 1 ,Y t - 1 ), 

C t -1 Et -1 

-~V-" 

predictor 

( 13 ) 

where (a) is obtained with the Markovian property and 
Bayes rules, and the current observations Y t and Ut work 
as correctors to improve the prediction of the current states 
based on the past. Inside ( fT3) , the probabilities of the corrector 
part could be computed by the emission matrices from the 
parameter vector rf , the first two probabilities of the predictor 
could be computed by the transition matrices in rf, and the 
last probability is the previous conditional probability of the 
hidden states, which is obtained from the previous step in 
our recursive algorithm. Note that the recursion procedure is 
exactly connected by the last probability term inside ( fl3| . as 
it is the only term dependent on the previous information. 

After we compute ( fl~3| >, to obtain the objective in GD. 
the value of P{Ut,Y^U t 1 , Y t 1 ) has to be calculated 
in order to divide it out from (fl3| . Thus the second 
step calculates P(U t , y t_1 jT - As the summation of 

P(Ct, EtlU*, Y*) over all possible hidden states equals one, 
the value of P(Ut, Y t_1 ) could be calculated as 

p(u t ,y t |u t - 1 ,y*- 1 ) 

= P(U t , YtfU*- 1 , Y*- 1 ) Y Y P ( Ct > E *\ ut . Y*) 

Ct E t 

= YY P ( Ut ’ YtlU^^-^PiCt, E t \U\ Y*). ( 14 ) 

Ct E t 

Although the individual probability terms inside (|T4} could 
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not be directly computed, the multiplication of the two proba¬ 
bilities actually is already the result we obtained in (jT3j in our 
first step. Then marginalizing over all possible ( C t ,E t ) pairs 
leads to the conditional probability of current observations 
P{Ut,Y t \U t ~ 1 as we desire. 

At the last step, the conditional probability of hidden states 
P(Ct, E^U*, y*) is calculated by dividing ( p~3] > with 
which actually could be treated as a normalization procedure. 
It is worth noting that when t = 1, which is the starting point 
of our recursion algorithm, the multiplication corresponding 
to is expressed as 

PiUuYjPiEuC^UuY!) 

= P(E 1 ,C 1 ,U 1 ,Y 1 ) 

= P(U l ,Y l \E 1 ,C 1 )P(E ll C 1 ) 

= P(U 1 \E 1 )P(Y 1 \E 1 ,C 1 )P(C 1 \E 1 )P(E 1 ), (15) 


where the first two probabilities are from the emission matrices 
in the parameter vector rf and the last two probabilities are 
from the initial probabilities n c and n E inside rj. Using the 
same method as in < p~4] ) to obtain the value of P(JJ\. Y )) and 
dividing the multiplication result in ( fl5] > by P(U\,Yi), the 
conditional probability of the hidden states at t = 1 could be 
computed. 

After obtaining the conditional probability of the hidden 
states P(Ct, EfU 1 ,Y t ), the maximization procedure in © 
could be solved by a simple 2-D exhaustive search, since the 
cardinality of the possible hidden states are within a finite 
range. 

In summary, the algorithm to compute the conditional 
probability P(Ct, EtlU*, Y*) is summarized below. 


Algorithm to Compute P(C t , E t \U l , Y l ) 
Input: ?y, U t , Y* 

1 for * £ .5? 

2 for j £ if 

QWK 


P(E 1 = i,C 1 =j\U 1 ,Y 1 ) = 


E E 

■ £S£ 3<=<iS 


4 end 

5 end 

6 for k = 2, • • • , t 

7 for i £ if 

8 for j £ E 

9 P(E k = i,C k =j\U k ,Y k ) = 

Q(k) E E bmj,iaiiP(Ck — i=m,Ek-i=l\U k - 1 ,Y k ~ 1 ) 

_ im _ 

E E Q( fc ) E E b m j, i ai i P(C k - 1 —m,E k _ 1 =l\U k - 1 ,Y k ~ 1 ) 
o=s? iz’e u=se 


10 end 

11 end 

12 end 


where Q(k) stands for dnj k - , if -. 

Remark: For an energy harvesting based PU, there is always 
a collection of energy states that are not enough to support 
reliable PU transmissions. Let 2z?o C 2z? denote the subset 
that contains all the energy states from ,E that are insufficient 
for reliable PU transmissions and Jzfj be the complement of 
..2)i. From the relationship given in <(TJ, when E t £ 2z?i, the 


estimated channel state C t could be either 1 or 0. However, 
when Et £ Ct could only be zero. This is due to the fact 
that the joint probability of C t = 1 and E t £ 2z? 0 is always 
equal to zero. As such, the optimal solution for could not 
be (Ct = 1 ,E t £ -2o)> which could be eliminated from the 
searching space in solving ( |TTj ) to reduce the computational 
complexity. 

B. Performance Comparison via Mutual Information 

Compared with our proposed 2-D sensing scheme, the tradi¬ 
tional PU detection methods do not take the spatial correlation 
of energy harvesting processes into consideration, i.e., for 
traditional methods the only available observation is Y f . It 
implies that under the same MAP principle, the traditional 
sensing scheme is basically solving 

C t , E t = argma xP(C t = c t ,E t = e t ,y 4 ). (16) 

c t ,e t 

The mutual information between the hidden states and 
observations could be used to quantify the estimation perfor¬ 
mance. In our setup, it can be evaluated as 

I(C t , E t -Y\ U f ) = I(C t , E t - F‘) + I(C U E t -U^Y*), (17) 

where I(Ct, E^Y 1 ) is the mutual information quantifying 
the performance of the traditional detection methods, and 
I(Ct, E t :, U* |Y 4 ) is the information gained by additionally 
considering the power domain information. The mutual in¬ 
formation gain is zero if and only if (C' t , E t ) are independent 
of U 4 conditioned on Y*, which is not the case in our system, 
such that the mutual information gain is always greater than 
zero in our method. 

The above mutual information gain demonstrates that U * 
could help the estimation of the PU hidden states; therefore, 
the proposed 2-D sensing method has the potential to perform 
better than the traditional detection methods, which will be 
further validated by simulation results in Section VI. The 
benefit of using energy observations for PU detection could be 
much more remarkable when the PU-to-SU links experience 
deep fading or shadowing, where the channel observations are 
much less reliable but energy observations are not affected by 
such channel impairment. Simulation results in Section VI will 
further illustrate this point. 

V. Parameter Estimation 

In this section, we discuss the proposed algorithm for pa¬ 
rameter estimation. We need to first define several intermediate 
variables, which are used for algorithm implementation, then 
derive the algorithm based on the traditional EM algorithm. 
Note that the proposed algorithm is off-line, with observations 
U T , Y 1 for training. 

A. Intermediate Variables 

For U* and E t , define 

a Y{t) = P(U\E t =i\ri), t = 1,2,3,- •• ,T, (18) 

which is the joint probability of getting the observation U* 
and having the PU energy state at time t ( E t ) equal i. 
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conditioned on 77 , with // defined in Section III. Similarly, with 
the communication channel output observation, define 

a[ j (t)±P(Y t ,C t = i,E t =j\r 1 ), t= 1,2,3, - - ,T, 

(19) 

which is the joint probability of getting the observations Y l 
and having the channel state C t as i and the energy level E t 
as j, conditioned on // as well. 

Since we estimate two hidden states with different observa¬ 
tions, the information carried by the joint observations should 
also be considered. Therefore, define 

(x u /{t)±P{U\Y\C t =i,E t =j\n), t = 1,2,3,- •• ,T, 

(20) 

as the intermediate variables for joint observations, which 
indicates the probability of getting the observations [/*, Y*, 
with Ct = i and E t = j. 

As both {E t } and {Ct} are Markov processes, with Markov 
property, the defined intermediate variable collections a u , a 5 
and o f ‘ V can be quantified recursively. For a u , the probability 
at time t, + 1 is given as 

af{t+l) = f{U t+ i\j)Y, a i( t ) a ij t = (21) 

iSJzf 

which is traced to <J{ (1) = nEiduji back at the initial 
state. This process is usually called a forward recursion, in 
which the next-time joint probability equals to the current joint 
probability weighted by the transition probability, followed by 
marginalization and multiplication with the emission probabil¬ 
ity. It can be interpreted that the forward recursion at each time 
equals to a predictor multiplied by a corrector. Similarly, the 
forward procedures for a Y and ) are respectively given as 

ajjil) = ncjinEjfiY^iJ), 

+ 1) = f(Y t +i\i,j) 'y } 'y } a m i{t)aijb m ij , 

t = 2,- ■ ■ ,T. (22) 

and 

alf (!) = ^Cji^EjdjuJ{Yi\i,j), 

0%;* (*+!) = djuJ(Yt+i\i,j) Y Y (I )<n j- 

t = 2, - ■ ■ ,T. (23) 

In addition to the forward recursion, backward recursions 
are then used to measure the probability of getting future 
output up to time T given the state at time t. Define 

py w = m +1 , u t+ 2, • • •, u T \E t = i, V ), 

PTj{t) = P(Y t+1 ,Y t+2 , ■■■ ,Y T \C t = i,E t =j,n), 

PV® = 

P(U t+1 ,Y t+1 ,U t+2 , Y t+2 , ■■■ ,U T , Y T \C t = i,E t = j , 77 ), 

t= 1,2,3, ••• ,T, (24) 

as the probabilities of getting the observation sequence from 
t, +1 to T, conditioned on the hidden states at time t. For the 


SU energy state, fj} ; (t) is calculated as 

PY(t) = Pi (* + 1 ) a ii d iu t+1 - * = 1, • • • , T - 1, (25) 
jesc 

with Pf (T) = 1, Vi. This indicates that the probability of 
getting the future observations up to time T conditioned on the 
current hidden states equals to such a conditional probability 
at the next time step multiplying the transition and emission 
probabilities of all possible current latent states. As it can be 
seen, the backward recursion uses the future observations to 
smooth the current inference. Meanwhile, the calculations for 
PYj(t) and pfj' 1 ( t ) could be expressed as 

PYj(t)= Yl YP^ t + l ^^ Yt C^ rn ^ a i lbim P 
Pij Y (t)= Y YP™i ( t + 1 )f( v t+i\mJ)di Ut+1 a j ib tm ^i, 

t = 1,- ■ ■ ,T — 1, (26) 

with p V(T) = 1 ,p% Y (T) = 1, Vi, j. 

After obtaining the expressions for forward and backward 
recursions, we need to define two update variables for each 
observation process to make the implementation of the pro¬ 
posed parameter learning algorithm more straightforward. One 
update variable is defined to consider the joint probability of 
getting the entire observations from time 1 to time T and 
the hidden state at time t\ the other one considers the joint 
probability of getting the entire observations and the hidden 
states at time t and t + 1. Specfically, for the SU energy 
observation process, define 

1 Y{t) = P{U T ,E t = i\ V )=aY{t)l3Y{t), t = 
eY j (t)±P(U T ,E t =i,E t+1 =j\r ] ) 

= aY {t)aijl3y [t + 1 )d jUt+1 , t = 1, • • • , T - 1. 

(27) 

Clearly, the update variables above combine the forward and 
backward variables with transition and emission probabilities. 
For the channel observation and the joint observation processs, 
variables are defined as 

7 Yj(t) = P(Y T ,C t = i,E t =j\v) = aYj(t)/3Yj{t), 

£ i7,mz(0 - p (. YT . c t =i,Et= j, C t +1 = m, E t+ 1 = l\rj) 

= a Jj(t) a jlbim,lPml(t + 1 )/(^t+l|TO, l). (28) 

and 

7 % Y (t) ± P{U T ,Y T ,C t = i, E t = j\rf) = 

Inlit) ~ P(.U T , Y t , C t = i,E t = j, C t+ 1 = to, E t+ 1 = l\rj) 

= a(t)a j ib imil pY;Y {t + l)f{Y t+1 \m,l)di Ut+1 . 

(29) 

Again, the reason to define these intermediate variables is 
that they will be used to implement the proposed parameter 
estimation algorithm for HIMM, which is introduced in the 
following subsection. 
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B. Parameter Estimation Algorithm 


The EM algorithm ifZTl is adopted here used to iteratively 
find the unknown parameters that maximize the likelihood in 
a statistical model, in which the latent variables exist. Here, 
E stands for expectation and M stands for maximization, 
and the algorithm alternates between the E and M proce¬ 
dures. The current expectation step creates an expected log- 
likelihood function, which is averaged over the latent variables 
given the estimated parameters in the previous step; then the 
maximization step is designed to find a new parameter that 
maximizes the function of log-likelihood formulated in the 
current expectation step. The algorithm continues until the 
results converge (when the difference between the results of 
two consecutive iterations is within an acceptable region). As 
proved before, the expected log-likelihood is non-decreasing 
during the iterations and the algorithm converges; however, 
it is not guaranteed to converge to a global optimum. This 
is due to the fact that although the maximization step is 
seeking the global optimum for the expected joint probability 
of hidden states and observations, the results may not be the 
global optimum for the likelihood function, which is from 
the conditional probability of observations. Detailed analysis 
about the EM algorithm convergence could be found in l22l . 

1) Expectation Step: In the expectation step, we need to 
consider the average of the log-likelihood function over the 
latent Markov states. If, in an extreme case, the Markov states 
are known, then the E step could be decomposed into T 
number of subproblems since the temporal dependence of 
the Markov chain contains no additional information for the 
expectation of the log-likelihood function. Let rjl k ~ 1 '> denote 
the parameter estimated at the previous iteration; and let 
O = {U T ,Y T } stand for the collection of observations. The 
expectation of log-likelihood can then be expressed as 


Hm v ik ~ 1] ) = E C T iBT {iog p(e t , c T , o\ v )} 


= E E P(E T ,C T \0, V ( k ~V) log P(E t ,C t ,0\ V ) 

E T(zcf T c T e c e T 


= E E 

E T(zjf T c T ev T 


P(E T ,C T ,0\r ]( k - lS >) 
P(0| ?7 (fe- 1 )) 


log P(E T ,C T ,0\ V ). 


(30) 


Inside the equation, E and E are abbreviations for 

E T ei? T C' r e < r r 

E E • • • E and E E • • • E ■ The achieved 
E 1 e^‘ -E 2 e-S? e t g. s? c 1 o c gc 2 ^ c T ev 

expected log-likelihood is utilized in the Maximization step to 

find a better parameter estimate, which should not decrease 
the averaged log-likelihood. 

2) Maximization Step: The objective in the M step is 
to achieve a better parameter estimate r )( k \ which is con¬ 
sisted of {n E ( k \ Tt c ( k \ A^ k \ D( fe ), n^ k \ <r 2 ( fe )}. Inside 
L(rj', r)( k ~P), since the denominator P(0\p^ k ~^) is a constant 
with respect to tj^ k \ the M step only needs to maximize 

L {n\ t7 (fe_1) ) / = 

E Y P(E T , C T , 0 |? 7 (fc — 1} ) log P(E T , C T , 0\rj) 

E t £& t c T <z c e T 

( 31 ) 


over the parameter. However, based on the structure of Fig. [T] 
the log-likelihood can be decomposed as 

logP(E T ,C T ,U T ,Y T \r)) 

= log {P(U T , Y T \E T , C T , V )P(E T , C T \r))} 

= \og{P(U T \E T , V )P(Y t \E t , C T 1 if)P(C T \E T , V )P(E T \r,)} 
= logP(U T \E T ,rj) + logP(Y T \E T ,C T ,r]) 

+ log P{C T | E t , p) -flog P(E t 1 77 ), (32) 

where each element depends on different optimization decision 
variables and there are no decision variables that have been 
shared by any two elements above. This implies that the 
maximization over the entire log-likelihood could be factorized 
into several independent subproblems, while the summation 
over all the optimal subproblem results leads to the optimal 
result of the entire problem. Mathematically, it indicates 

maxL(r?; 

= max E E p (E T ,C T , 0 \r 1 {k - 1) )lo g P(U T \E T ,r 1 ) 

E T e cffT C r g¥ T 

+ ma X E E P(E T ,C T , 0 \ V (k ~ 1) )logP(E T \r 1 ) 

n ' A ET eJSfT C T 

+ max Y, E p (E T ,C T ,0\ V ^ k - 1) ) log P(C t \E t , V ) 

” 3 ET C T 6¥ T 

+ max Y E P (E T , C T , 0\p (k ~ 1) ) log P{Y T \E T , C T , p). 

IJ " rT ~ E T eZ£ T c t gv t 

(33) 

By solving the four subproblems above, the optimal solution 
for could be achieved, which is then used for the next 
iteration if the convergence requirement is not yet satisfied. 

a) Optimal D :/,: L First, let us focus on learning the 
parameter D, which is the emission probability of the PU 
energy arrival process. Since D = {d^}, learning d, :) for 
all i, j £ .Y J is sufficient to have the estimated value of D. 
Before moving on to the maximization procedure, the objective 
function could be further simplified as 

max Y P (E T , U T \’q ( ' k ~ p ) log P(U T \E T , 77 ). (34) 

E T e£f T 

This is due to the fact that the hidden channel states C T could 
be marginalized by summing over all the possible outcome, 
and also 

P(E t ,U t ,Y t \ V ^) 

= P(U T , e t \y t , 7? ( fe - 1 ))p(y T | 7? ( fe - 1 )) 

= P(U T \E T ,rt k ~ 1) )P(E T \Y T ,p (k - 1) )P(Y T \ii k - 1) ), 

where P{Y T \r) ( ‘ k ~ 1 '>), P{E T \Y T ,^ k ~ p ) are constant values 
with respect to D. Besides, P(E T \r) < ' k ~ 1 ' > ) is positive and inde¬ 
pendent of the decision variable D, such that its multiplication 
with the objective function will not affect the estimation of 
D. As such we consider joint probability P(E T , U T \r] ( ' k ~ 1 ' > ) 
instead of the conditional probability P(U T \E T , rj^ k ~ 1 '>) with¬ 
out influencing the maximization step. 

On the other hand, since dij is the emission probability, 
the objective function needs to be decomposed from a joint 



probability into T independent emission probabilities, in order 
to estimate the value of dij. The decomposition procedure is 


E P(E T , U T \r^ k ~ r> ) log P(U T \E T , rf) 

ET<zS£ T 

= E 

E T e& T t=1 

T L — l 

= EE E P(U T , E T \r 1 (k - 1) )5{E t - i) \ogP(Ut\Et = i, rj) 

t=1 i =0 E T e& T 
T L — l 

= E E p ( [/r - E < = *l^ _1) ) logP(U t \E t = i, V ) 

t=1 i=0 

T L—l L—l 

= E E P ( [/T > Et = *l’l ( ' C_1) ) lQ g E 5 <P* - fidij- (35) 

t=1 i =o j —o 

Inside the equation, 5(x) is the indicator function such that 
8(x) = 1 when x = 0 and 8(x) = 0 for any other values of x. 
Since dij denotes the emission probability of the PU energy 
level, for any i the summation of d, t] over all the possible 

L -1 

outputs equals to one, i.e., dij = 1, Vi. Therefore, the 

j—0 

maximization problem is finalized as 

T L — l L—l 

m D X EI El = *|?? (fc_1) )l0gE S ( U t ~j)dij 

t= 1 2—0 J —0 

L—l 

s.t. E ^j = 1 (36) 

j=o 

This is a typical optimization problem to find a maximum of 
the objective function subject to an equality constraint, which 
could be solved by the Lagrange multiplier method lf23l . After 
calculations, the optimal result for d\j is given as 

J2 6(U t -j)P(U T ,E t = i\r,( k -V) 

k ) _ L=1 _ 

jrP(U T ,E t = i\r 1 ^-P) 

t= 1 

E 8{Ut-jh¥(t) 

= -^-, i.jGJSf. (37) 

Et Y(t) 

t=l 

which implies that the estimated probability of emitting j at 
state i for the PU energy level equals the number of times 
that the output from the latent state i is j divided by the total 
number of times that the latent state is i. 

b) Optimal n E ( k \ A^ k \- The objective function for this 
subproblem could be further decomposed into two parts, 
in which the decision variables are n E and A = {a t J }, 
respectively. The decomposed objective function is 

E E p (E T ,C T ,U T ,Y T \i 1 ^)\ogP(E T \r ] ) 
B T ei? T c T e^ T 

= E P(E T ,U T ,Y T \^ k ~ 1 '>)\ogP{E T \i 1 ) 

E T e^ T 


= E P(E T ,U T ,Y T \^ k -^) 

e t gs? t 

x jlogTT^ +E l0 S a -Et-lSt| ■ ( 38 ) 

Clearly, the optimization procedures to maximize the objective 
function over n E and A could be treated as two independent 
processes. Therefore, after factorizing the joint probablities 
into T independent probabilities based on the Markovian 
property, the problems become 

max V P([/ T , Y t , E l = i\p {k ~ 1] ) log nf 

7T E ^ - J 

L -1 

s.t. E Ttf = 1) i G Jzf, (39) 


and 


:a. 


tEEE P{U T ,Y T ,E t - 1 =i,E t =j|r ? ( fc - 1 ))l 0g 

a%D t =2 itS£ j 
L -1 

s.t. E a ij — * G . (40) 

3=0 

By solving the subproblems with the Lagrange multipler 
method, the estimated values of tt e and A for step k are 

^(fc) = P(U T ,Y T ,E 1 =i\ V < k - 1 '>) 


2 J 


P(U T ,Y T |r/( fe-1 )) 


e i-r (1) 


jev 


and 


(k) t =2 

a ij — T 


E E Hj (!) 

jetfieS? 


E P{U T ,Y T ,E t - 1 = i,E t =j\rt k ~V) 


(41) 


E P{U T ,Y T ,E t _ 1 =i\i 1 ( k - 1 '>) 

t =2 


T 

E 


E E 


E E 

t =2 me? 


(42) 


for any i £ Jzf. It can be seen that as channel and energy 
observations contain the information of the PU energy level, 
they both contribute to the estimation of the PU energy 
transition matrix. 

c) Optimal Tc c ( k \ B( fc ).- Following the previous argu¬ 
ments, the objective function here can also be decomposed 
as 

E E P(U T ,Y T ,C T ,E T \ v ik - 1 ) )\ogP(C T \E T ,r 1 ) 

E T eS g T C T e ^T 

= E E P{U T \E T ,^ k ~ 1 ) )P{E T ,C T ,Y T \r 1 P- 1) ) 


ETeseT c T ev T 


log7T G +E 1 °g b C? t _ 1 Si,C t >• 


(43) 
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Since P{U T \E T ,? 7 ^ fc_1 ' ) ) is not changing with respect to tv g 
and B, it can be neglected for estimating these parameters. 
After factorizing the objective function into T independent 
elements, the optimal parameters are obtained as 


^ k) = 


P(Y T ,E 1 =i,C 1 =j\r ,( k - 1 )) 


P(Y T ,E 1 =i\ v ( k ~ 1 '>) 


t£(i) 

Et£(i) ! 

jev 


i £ SY ,j £ c to , 


and 


(44) 


E p (Y T , Ct- 1 =i,E t = k, Ct = j\v {k - 1] ) 

b (.k) = t =2 _ 

f;P(Y T ,Ct-i=i,Et = k\ri<. k - 1 )) 

t= 2 

E E eL—!) 

= ;- 2 -T° -. (45) 

E E E eLjfcft- 1 ) 

£=2 m =0 


c/) Optimal pi, a 2 : The channel output {V)} is a conti- 
nous states process, such that given the latent states of channel 
and PU energy, { Y t } follows a Gaussian distribution for any 
t. Learning p = {ptj} and variance <r 2 = {of-} for i £ ‘if 
and j £ if of the conditional Gaussian distribution could 
determine the emission probabilities of channel and PU energy 
states. 

Conditioned on the hidden PU energy process, the obser¬ 
vation of the SU energy level is independent of the hidden 
channel state and the channel output; therefore, U T does not 
contribute to the estimation of the channel output parameters. 
Since given the latent states, the channel output {Y t } is 
independent over time, we have 

E E P(Y T , c t , e t |» 7 <fc - 1 >) io g p(y T |£ T , c T , n) 

E t C t 

= E E P ( yT - ° T > E T \v (k - 1) ) log fl P(Y t \E t , Ct, r,) 

e t c t ‘=i 

= E E E p ( yT ’= <. ■ e t= jW k ~ 1] ) 

t=i iev je £f 

x log P{Yt\Ct = i, E t = j, f}), (46) 

where P(Y t \C t =i,E t = j, r]) ~ N(Y t ; Pij,a^) and E- E 

e t c T 

stand for E > E respectively. Since this is an un- 

e t gj? t c T e c e T 

constrained convex optimization problem with a differentiable 
objective function, letting the first derivative of the objective 
function equals to zero and solving the rest will provide the 
optimal point for pij and of-. Then the estimated values of 
these parameters are 

E P(Y T ,C t = i,E t = j\^ k -P)Y t E 7 TjWt 

(k) _ t=l _ _ t=l _ 

E P{ Y T 1 Ct =i,E t = j\ V ( k ~P) E 7 lit) 

t= 1 t—1 

( 47 ) 


and 

E P{Y T , Ct = i,E t = j\^ k ~P)(Y t - p^) 2 

2 (fc) t= 1 

(j.) = - 

J2P(Y T ,C t = i,E t =j\r 1 ( k -V) 

t= 1 

Et 

= —- T -■ (48) 

Et lit) 

t—1 

However, the expressions above only work for i = 1, Mj £ ,Y. 
Since when i = 0, PU does not transmit signal and the channel 
output does not contain any information about the PU energy 
level. As we mentioned, in this case /j i; and variance of- are 
identical for any j £ ,Y. Then the estimated values of these 
parameters are 


and 


(fc) _ t —l 
Pij T 


E P(Y T ,Ct = i\v {k - 1) )Y t 


J2P(Y T ,C t =i\r ] (fc-D) 

t= 1 


T L — l 


E Et m>Y t 


t =1 I j =0 


T L—l 

E E 7&(t) 

t =1 J =0 


(49) 


E P{Y T ,C t = i\^ k ~P){Y t - ^ fc) ) 2 

CT 2(fc) = t=i__ 

j:P(YT,C t = i\r 1 ( k -P) 

t=l 

_ t=l [j =0 J 

— T L—l ’ 

E E 7^(0 

t=i i=o 


(50) 


when i = 0 . 

Based on the calculations above, the estimated parameter 
could be updated to ij lk> . As the algorithm is derived from EM 
algorithm, the expected log-likelihood is non-decreasing over 
iterations. Let e be the acceptable difference for convergence, 
if L(rp,r ) (fc) ) - L{rp, ? 7 (/c_1) ) < e, we claim ?/ fe ) as the fj na l 
estimated value for parameter tj and use it for further appli¬ 
cations; otherwise, we use r/Y to calculate the expectation of 
the log-likelihood and continue EM iteration to get ij <k+ i ’. 


C. Initialization 

As can be seen from l24l . since at every iteration the 
estimated parameter increases the expected log-likelihood 
function value, the proposed algorithm based on the EM 
algorithm converges to a stationary point, which can be a 
local optimum or a saddle point. In other words, the global 
maximum is not guaranteed such that convergence to which 
local optimum depends on the initialization of the algorithm. 
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x iq 4 log likelihood 




Fig. 3. Increasing log-likelihood with Parameter Learning pig. 4. Comparison of Detection Performance 


For simplicity, a multi-try random initialization routine is 
adopted in our algorithm. Since at each time the algorithm 
with a random initialization converges to a local optimum, it 
is logical to run the algorithm multiple times with different 
initializations and choose the outcome that has the maximum 
likelihood value. Thorough discussions of optimal initializa¬ 
tion strategies is out of the scope of this paper. 

In summary, the proposed algorithm for the HIMM param¬ 
eter estimation is illustrated below. 


Parameter Estimation Algorithm 


1 

2 

3 

4 

5 

6 

7 


Define iteration variable k = 1. Take the initialization 
value t/ 0 ) = {n E ^, A(°>, B<°), D(°>, ji (0) , <r 2(0) } 

from certain distributions. 

Use to calculate intermediate variables 

a u ,p u ,f,e u ,a Y ,p Y ,i Y ,e Y ,a u ’ Y ,/3 u ’ Y , 1 u ’ Y ,e u ’ Y . 
Use the intermediate variables from step 2 to calculate 
the current estimated parameter rj^ k \ 

Based on calculate the expected log-likelihood 

L{j]\ t? (fc) ). 

Repeat steps 2 to 4 until L{r]\ r — L(r ]; r < e. 
Repeat steps 1 to 5 with different initializations to find 
r/' 1 ' 1 with the largest likelihood value. 

The estimated parameter is set rj = 


VI. Simulation Results 

In this section, numerical results are provided to validate 
the effectiveness of the proposed parameter learning and 2-D 
sensing algorithms. Our scheme considers both channel and 
energy observations while the reference energy detector based 
spectrum sensor, which is one of the well-accepted traditional 
PU detection methods, only uses channel observations. 

For the parameter learning, we consider a training model 
with 5000 channel outputs and SU energy observations. The 
trend of the increasing log-likelihood with the proposed pa¬ 
rameter learning algorithm is showed in Fig. [3] It can be seen 
that there is a sharp log-likelihood increase in the first five 
iterations, and after 60 steps the algorithm converges. As one 


property of the EM algorithm, the log-likelihood converges 
monotonically and every random initial value guarantees the 
convergence. In this simulation, we tried 15 different initial 
values, and the estimated parameters are obtained by choosing 
the one with the maximum converged log-likelihood value. 
The initial values are randomly chosen according to the uni¬ 
form distribution between 0 and 1 with proper normalizations. 

In Fig. [4] we show the differences of probability of detection 
versus SNR at the SU between the traditional sensing and 
the proposed 2-D sensing methods, with HIMM parameters 
learned by the proposed learning algorithm. The SNR at the 
SU is defined as the SU received signal power (from PU) 
divided by the channel noise power. In this figure, the channel 
noise level is set as = 3 dBw, and the channel path loss 
is —4 dB. The energy state takes values from {1,2,3,4}, 
in which the insufficient energy state subset is Jz?o = {1}. 
Note that with each SNR value, we set the two schemes in 
comparison to have the same false alarm performance. 

From Fig. [4] we see that the 2-D sensing method outper¬ 
forms the traditional energy detector based scheme over all 
SNR values, since the 2-D method considers observations 
of both the channel state and the PU energy state, and it 
also utilizes the hidden Markov signaling structure between 
the PU and the SU. As the estimator used to estimate the 
hidden channel and energy states is a MAP estimator, it 
actually provides the optimal result that jointly estimates the 
hidden states of the Markov model. When the SNR is low, 
the observation over the PU-to-SU channel is not reliable for 
PU detection, which largely affects the performance of the 
traditional method. Since 2-D sensing also takes the energy 
observation into account, which is not affected by the PU-to- 
SU channel, the performance of 2-D sensing is much better. As 
SNR goes up, the channel observations become more reliable 
and the benefit of using the additional energy information is 
less, which implies that the relative gain of using 2-D sensing 
will decrease. However, as the 2-D sensing method does 
not sacrifice any information during the estimation, although 
the gain may become less, 2-D sensing still outperforms the 
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Fig. 5. Tracking Performance of Energy States 


traditional method in the high SNR region. 

Another advantage of using the proposed 2-D sensing 
method is that SU could estimate the transmit power of 
PU, which is proportional to E t for broader CR applications 
mm. In Fig. [5] we show the tracking performance of the 
hidden PU energy state E t , where it can be seen that the 
proposed method could estimate the hidden energy state quite 
well. 


VII. Conclusion 

In this paper, we applied a HIMM model to represent 
the interaction between energy harvesting based primary and 
cognitive radios, and provided an EM based algorithm to 
estimate the parameters in this HIMM structure. Then we 
proposed a novel 2-D sensing scheme, which jointly considers 
the observations from the spectrum and power dimensions. 
The proposed scheme could sense the spectrum and estimate 
the energy level for PU transmission simultaneously. We 
showed that the proposed 2-D sensing method outperforms 
the traditional spectrum sensing method, since it utilizes the 
facts that the PU power usage and channel usage are inter¬ 
dependent events, and the PU/SU energy harvesting processes 
are spatially correlated. 
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